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We study numerically the decay of the inflaton by solving the full non-linear equations of motion 
on the lattice. We confirm that parametric resonance is effective in transferring energy from the 
inflaton to a scalar field as long as the self-interactions of the second field are very small. However, 
in the very broad resonance case (q 1) the decay rate is limited by scatterings, which significantly 
slows down the decay. We also find that the inflaton cannot decay via parametric resonance into a 
scalar field with moderate self-interactions. This means that the preheating stage may be completely 
absent in many natural inflationary models. 
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I. INTRODUCTION 



According to inflationary scenarios of cosmic evolution the Universe undergoes a period of exponential expansion 
after which it is essentially devoid of matter. At this stage the energy density is almost entirely in the large oscillating 
expectation value (EV) of the scalar field that drives inflation, the inflaton. The theory of reheating is concerned 
with the question how this energy is transferred from the inflaton to other fields and eventually thermalized. Early 
calculations assumed that the inflaton would decay to lighter particles via perturbative decay ffl . It was only recently 
realized J^] that non-perturbative mechanisms may completely dominate the reheating process. In particular the 
classical phenomenon of parametric resonance may lead to explosive particle production and very rapid decay of the 
inflaton field ||]. This stage of the evolution, called preheating in the literature, leads to very different physics than 
would be obtained from slow perturbative decay. Different aspects have been investigated by many authors [Q . 
In this paper we present the first fully non-linear study of inflaton decay into a second scalar field. By fully non- 
, linear we mean that we take into account the back reaction of the created particles on the equations of motion as well 
as their scattering. This is accomplished by integrating the classical equations of motion of the system with certain 
random initial conditions, a technique that was first applied to the one field case in Q . The system is modeled by a 
■ simple two scalar theory with effective potential 
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We will focus on two general classes of inflaton potentials V((f>): type I for which (f> <C mj \/Xp (for these we will set 
\<p = 0), and type II for which cf> ^> m/VX^ (for these we will set m — 0). In both cases we will assume that there is 
a channel for the inflaton to decay into a light scalar and thus set the tree level mass /i of the x-field to zero. The 
' effects of a massive decay product have been investigated using the Hartree approximation in |6| . For type I models 
we will neglect the expansion of the Universe M . For type II models our calculations are directly applicable to the 



expanding Universe, as will be explained in section IV 



With the above approximations and neglecting A x the equations of motion for the modes of the x-field are 



d 2 Xk 



+ (k"+g4>"(t))xk=0. (2) 



dt 2 

For type I models, where the oscillations of the </>-field EV are sinusoidal, this can be written as the Mathieu equation 

d 2 Xk 



dz 2 



+ (A{k)-2qcos(2z)) X k = Q, (3) 



where z ~ oj^t, A(k) = k 2 /cj? + 2q, q = g$ 2 /4lu 2 , lo^ = m is the frequency of the oscillations, and $ is the slowly 
varying amplitude of the inflaton EV. For type II models the oscillations of 4> are given by elliptic functions, but for 
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the illustrative purpose at hand it is sufficient to replace the periodic EV by a sinusoid of the same frequency. In that 
case the \ modes again satisfy (||), except that now uj^ = c^/X^, where c ~ 0.85. The important point is that @ has 

exponentially growing solutions of the form \k oc exp(/ii 2) within a set of resonance bands labeled by the integer 

index n. This growth of the modes corresponds to exponentially growing occupation numbers oc exp(2/x^ u)st) 
and may be interpreted as particle creation. Parenthetically we remark that for type II models the inflaton may also 
decay into its own fluctuations. Writing <f> = (cf>) +5<j> and linearizing one finds that the S<pk also (approximately) obey 
(|J), with A w fc 2 /cj0 + 2.08 and q sw 1.04. For g ^> A^ the decay into ^-fluctuations dominates, however, because the 

relevant values of /z£ are a few times larger. 

The literature distinguishes between narrow resonance, defined by q <C 1, and broad resonance, for which q > 1. 
While the first case can be analyzed analytically Q], it is more difficult to obtain quantitative results for the broad 
resonance case ||. It is precisely this regime, however, which is most interesting: models generally have q ^> 1 at the 
end of inflation, and it is broad resonance that leads to explosive particle creation. To get a feel for the values of the 
parameters involved, note that typically <!>o ~ -Mp at the end of inflation, where Mp = 1/V87tG = 2.4 x 10 18 GeV is 
the reduced Planck mass. COBE data || restrict m < 10 -6 Mp for type I models and A^ < 10~ 12 for type II models. 
To prevent radiative corrections from generating an inflaton self-coupling in conflict with COBE, one generally needs 
9 2 ^5 ^(t>\max, i-e. g ;$ 10 -6 . Finally, if x is to represent a "standard" field then there is no reason for its self-interaction 
(as well as its coupling to other fields) to be tiny. Thus one might expect A x ~ 10~ 2 — 1. The upshot of all this is 
that there is a "natural" hierarchy for the couplings in the model, 

A x > g > A , (4) 

and that q is generally large at the end of inflation. We believe that the role of A x is of particular importance, especi ally 
since the final state self-coupling has been ignored in much of the literature so far. The exception is reference [p.0| , 
where it is was found that the self-interactions can have important effects in the case of narrow resonance. Below we 
will see that the same holds true in the physically important broad resonance regime. 

The difficulty in analyzing the broad resonance case stems from the fact that all parameters in the Mathieu equation 
for the modes vary quite rapidly. For q very large particle production takes place during a tiny fraction of the period 
as the EV passes through zero. The amplitude $ of the oscillations decreases as energy is transferred from the zero 
mode to the fluctuations. The produced \ particles generate a contribution to the mass of the inflaton of the form 
g{x 2 )t which changes to^. The x-field self-interaction, which we ignored in deriving the Mathieu equation for the 
fluctuations, produces an effective \ mass of the form m 2 eff « 3A x (x 2 )- This adds to A{k) the time dependent term 
m 2 gg/w 2 ,. In addition to these backreaction effects there is scattering: the resonance produces particles in narrow 
momentum bins with large occupation numbers, and one certainly expects these to scatter and spread out rapidly. 
This effect turns out to be very important, and it has not previously been studied in two field models. 



II. THE METHOD 



In order to take all of the above into account quantitatively we discretized the exact classical equations of motion 
and solved them on a three dimensional lattice. The classical equations of motion are a good approximation to the 
dynamics provided the mode amplitudes are large in the sense that the commutator of the canonical variables can be 
replaced by the Poisson brackets, which reduces to (f>k<t>k S> 1 |lfyj]. Although our initial state, which will be described 
below, does not satisfy this condition, it is satisfied very soon after the resonant decay begins. It is useful for the 
numerics to work in terms of dimcnsionless variables. We use x M = ^/g^^x^, rh = m/^/g$>o, A^ = X$/g, X x = X x /g, 

4> = 4>/$>o, and X = x/^o- In terms of these variables the initial inflaton EV amplitude $o as well as g completely 
disappear from the equations of motion. While $o simply sets the scale in the problem and becomes irrelevant for 
the dynamics, the coupling g is still important: it appears in the equations for the random initial fluctuations and 
regulates their amplitude compared to the EV. The fields are evolved using an explicit algorithm which is second 
order accurate in time and fourth order accurate in space. We varied both the number of grid points and the lattice 
spacing to check that we are near the continuum limit. The data presented in this paper were obtained on 128 3 
lattices. Regarding numerical accuracy we note that the total energy of the two field system was conserved to better 
than 1% in all of our runs. 

The initial conditions for the numerical simulation were chosen according to the following reasoning: the initial 
amplitude of the inflaton field is naturally of order the reduced Planck mass, $o ~ Mp = 1/\8ttG = 2.4 x 10 18 GeV 
[|12j . Typical values of the resonant momenta are given by A — 2q ~ y/q, which leads to 

(ss^y ^v^y-^^i^. (5) 
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For g > 1 this implies that k res /H ^> 1 at the end of inflation. Since this ratio does not decrease at later times, 
the resonant momenta are always on subhorizon scales. These momenta are not squeezed by inflation, and to a good 
approximation they are zero point quantum fluctuations. The initial conditions may then be approximated by a set of 
harmonic oscillators in the ground state. We define the phase space in terms of mode amplitudes {xk, Xk} for which 
(Xk) = (Xfc) = 0. These are defined as the minimum uncertainty states with (|xfc| 2 ) = l/(2u>k) and (|xfc| 2 ) = Wfc/2, 
where tUk is the appropriate mode frequency. In the semiclassical approximation these states are populated according 
to the phase space probability distribution. In our code we randomly generate gaussian distributed states with 
the width given by (2wfc) -1 / 2 for Xk and for Xk, and similarly for the <p modes. We point out that the 

evolution does not depend on details of the initial conditions, as long as the initial field configurations have the correct 
amplitudes p3| . 



III. MASSIVE INFLATON (TYPE I MODELS) 
A. Models with A x = 

Let us begin by discussing the simplest possible situation, namely a type I model with A x = 0. For the run presented 
here the parameters were chosen as follows: g = 1CP 8 , rn 2 = lO -12 ^, and lattice spacing Ax = bir/A. The initial 
amplitude of the 0- field EV, $o, sets the scale and does not effect the dynamics in any way, as discussed above. Note 
that for these values the q parameter in the Mathieu equation is q = 2500, putting us in the broad resonance regime. 

Figure 1 shows the EV of the inflaton field <j>, which starts decaying significantly around t — 5000 / ^/g^o = 0.5 x 
10 8 $o • Figure 2 shows the corresponding variances ((<5</>) 2 ) = {<ft 2 ) — (<j)) 2 and (x 2 ) of the fields. The variances increase 
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FIG. 1. The EV of the inflaton as a function of time (m = 10~ 6 $ , A x = A = 0, g = 10~ 8 .) 

roughly exponentially at first. This trend stops when 8m 2 , = g((5(f>) 2 ) becomes of order the resonant momentum 
squared. The reason is that the A(k) term in the Mathieu equation gets a contribution of the form g((8(f>) 2 ) /uii. 
When this term becomes of order fc 2 es /w| JT4J , the resonance band closest to A = 2q, which has the largest value of 
/ifc, becomes inactive. In the next band /ifc is down by a factor 3 |l5|] so the decay is slowed down dramatically. The 
rapid production of inflaton fluctuations, which kills the dominant resonance, may be surprising since there is no <p 
resonance in this model. These fluctuations are produced by scatterings of resonant x-Auctuations off the inflaton EV, 
which are fast, since both the <j) zero mode and x resonant mode amplitudes are large. It is interesting to note that 
if one analyzes this model using a Hartree type approximation which neglects scattering but includes backreaction 
effects, then the condition for the end of explosive particle production is that q(t) = g${t) 2 /4(w5+(7(x 2 )) becomes <C 1. 
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We emphasize that this is not the correct criterion: our simulations show that scatterings terminate the exponential 
stage long before q(t) < 1. For example, one sees from figures 2 and 3 that q(t) ~ 600 even at t — 5 x 10 8 $ ( J" 1 . 

The energy densities of the fields, as well as the total energy density of the system, is shown in figure 3. In 
the range t S [0.8,5] x lO 8 ^ 1 the field energy grows like t a , with a w 0.95 and a decay time scale of order 
Tdecay ~ lO 9 ^^ 1 . Based just on the Mathieu equation one would expect in this range exponential decay with time 
scale Tdocay ~ 1/ (2{fik)w<t>) ~ 10 7< i>Q 1 . The reason is that as <J> decays, the resonances sweep through the momentum 
space, and as they move to higher momenta, diminishes, and new ones emerge in the infrared (IR). In estimating 
the above time scale, we have taken the typical value for ~ 0.05. 

To understand why the actual time scale is much longer it is useful to look at the spectrum of fluctuations in 
k-space. To this end we define "occupation numbers" for the \~ field via 

uj x 1 

n k = -fiXkX.%) + ^pciXkX-k) > (6) 

where ui^ — \Jk 2 + g(4> 2 ), and the brackets denote averaging over directions. The nf. are defined similarly, except 
that ujf. = \/k 2 + m 2 + g(x 2 )- Figure 4 shows the occupation numbers of the x-field at various times. The growth 
is exponential at first, with the lowest three resonance bands clearly visible at k « .025, k « .09, and k « .13, 
respectively. The values of p^ obtained by assuming n£ oc exp u^t agree well with the predictions from the 
Mathieu equation. For the last time shown in figure 4, scatterings are already becoming important, and the resonance 
peaks get smeared out. 

The subsequent evolution is shown in figure 5. Scatterings are very fast after t ~ 0.5 x 10 8 $q , and the resonant peak 
structure is quickly washed out. The occupation numbers grow to about n\ ~ l/<?, and then remain approximately 
constant. There is a simple feedback mechanism which explains this behavior. Resonant particle production increases 
n£ , which increases the scatterings off the zero mode, which increases infrared 4> occupation numbers. This increases 
5m 2 , = g((S(f>) 2 ), which decreases ^ (cf. |l5[| ), slowing down resonant particle production and giving scatterings time to 
move particles toward higher momenta. This reduces ((5<fi) 2 }, which increases /x^, increasing resonant production and 
completing the feedback mechanism. Consequently n£ reaches an equilibrium value for which on average scatterings 
remove particles from the resonance as quickly as they are created, implying p^ os /p^ es — 2(j,kW<j> — r sca tt ~ 0. The result 
is a slowly varying distribution with a characteristic shape, as seen at late times in figure 5. The slope —kd\nn*/dk is 
small (< 2) in the infrared, grows as k increases, and eventually becomes greater than 4 at the scale A at which most 
of the x energy is concentrated. Particles created by the resonance diffuse to this scale A via scatterings and slowly 
move it to higher momentum. (Note that soon after scatterings become fast, A 3> k Tes .) Since the IR occupation 
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FIG. 3. The energy densities of the fields as a function of time (m = 10 _6 <I , o, A x = = 0, g = 10" 



numbers of the fields remain roughly constant during this process, the energy required is supplied by the zero mode. 
The number of <f> zero mode "particles" eaten up in moving one \ particle from k lcs to A is about Nq ~ A/u),/, 1, 
implying that scatterings are responsible for most of the EV decay. The upshot is that energy is drained from the 
zero mode at a slowly varying rate given by p'$ ~ Af^p^ 2/j,k^4>, where p* os is the energy density in the \ resonance. 
This implies a decay time scale 
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To obtain the last expression in Eq. (Q) we used 



d 3 k 



(2tt 



k 3 

—w n x 



where the resonance width w rcs is approximately given by 



0.C 



exp[-(m 2 +fc 2 )/(Afc) 2 



(7) 



(8) 



(9) 



as well as fc 2 0S — v 7 ?^ = v / 3 m </> (1) / 4 ' anc ^ ^ k — 0.23exp[-(m 2 + fc 2 )/(Afc) 2 ]. Here Ak ~ {y/gm^/2) 1 / 2 denotes the 
distance between the first two resonances. The numerical coefficient is obtained by noting that resonant energy produc- 



tion peaks around k ~ k rcs , and 



aim 2 



k 2 cs . The last relation remains reasonably well satisfied throughout 



the scattering regime. This is not surprising, since from 2pk^4> ~ r scat t one finds m 2 ~ fc 2 es ln(r scat t/2/zoW0), i.e. 
mr, is only logarithmically dependent on slowly varying quantities (from above, p,Q — 0.2/e). The quantities entering 
Eq. (Q) vary slowly during the scattering dominated regime, and one should take their typical values to get an order 
of magnitude estimate for the decay time. From figure 5 one reads off gn% ~ 0.04 and A/lu^ ~ 70, which leads 
to r d , 



ccay 



10 9 $q" . This agrees with the decay time scale seen in figures 1 and 3. The time scale in (|7|) is naively 
parametrically larger by ~ l/[ghi(l/g)] in comparison to the scale of parametric resonant decay. The reason why the 
decay time is not huge even for g very small is that grt^ ~ 0.04, i.e., the occupation numbers are ver y larg e. This 



however is not true when there is a moderate self-coupling of the x field, as will be discussed in section III B 



From the discussion above it is clear that the inflaton decay consists of two distinct regimes: fast exponential 
decay during which scatterings are irrelevant followed by a much slower decay governed by scattering processes. One 
can clearly separate these two regimes in figure 2: the scattering regime sets in at the time when the variances 
become slowly varying, which occurs at t ~ 0.8 x 10 8 $q Clearly there is a brief transition stage between the two 
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regimes, and since the exponential decay is a stimulated process, this is where most of the energy loss associated 
with the exponential regime occurs. The resonant occupation numbers stay roughly constant during this transition 
stage since the scatterings are already fast, but the infrared and near ultraviolet states are quickly filled. One can 
estimate quite generally what fraction of the inflaton energy has decayed at the time when the slowly varying state 
sets in, as follows. We have run our code for several values of g in the range g = 10~ 12 to 10~ 8 and find that 
g{x 2 ) — 9 J[d 3 k/(2-K) 3 u;^:}n^, ~ when the scattering regime begins This value is somewhat dependent 

on the exact initial position of the dominant resonance. Since at this time the variance is dominated by modes 
with k 2 ~ y/gui^o, we can estimate the occupation numbers of the \ field around the resonant momentum to be 
grv^ ~ 5^0/ y/g^o- At this time the x energy is dominated by approximately the same scale (y/guj^o) 1 ^ 2 , implying 
Px/(P4>)o ~ ^ip/y^^o = (l/2)?^ 1 ^ 2 . This means that for q 1 only a small fraction of the inflaton energy decays 
during the exponential regime fllTf , 

After the scatterings become important, the distribution broadens significantly, but the occupation number at 
the resonant scale remains to a good approximation constant. With our above estimate for gir£ we can thus 
give a parametric expression for the decay time in the scattering regime based on Eq. (Q). Our numerical results 
indicate that occupation numbers drop approximately exponentially between fc rcs and A (cf. figure 5). We take 
n k — n k ex P [ — (k — feres )/«]) where 1/k is the slope of lnrt£. As can be seen in figure 5, this value decreases slowly 
as the field decays. The relevant value of k can be computed by noting that when the EV decays significantly, about 
half of the original energy density (p^)o = ^$o/2 m hi the x fluctuations. This gives n 2 ~ q 1 / 4 v fg®ou<p/2. The scale 
A can be evaluated as A = (k) ~ 4k. Combining these results with Eq. (fjj) we arrive at the following parametric 
estimate for the decay time (recall that this is sensible only for q 3> 1 , since otherwise most of the field decays during 
the exponential regime): 

Tdccay ~3x 10 V — , <7 > 1 (10) 

The main point of this equation is that the decay time is very weakly dependent on q. (This conclusion is independent 
of our assumption that n£ falls off exponentially in the relevant momentum range, even though the exact power of q is 
not.) Note that surprisingly, for fixed and $o> the decay time actually slowly increases as g increases. That Tdecay is 
nearly independent of q in the scattering regime agrees well with our numerical simulations for 2.5 x 10 3 > q > 25. We 
have not verified this dependence for much larger values of q because such simulations require enormous computing 
times. We showed above that the fraction of energy that decays in the exponential regime is ~ l/V? oc Vv^' i m ply m g 
that paradoxically the field decays faster for smaller values of q oc g. Our numerical calculations indeed indicate that 
the fastest inflaton decay occurs for q ~ 1, rather than forg> 1. One should note, however, that in an expanding 
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FIG. 5. x-field occupation numbers vs. momentum at later times, when scatterings dominate (t* = -JgQot). 



universe parametric resonance is ineffective for q ~ 1 in type I models since the effective qt becomes much less that 
unity before the resonant occupation numbers grow large nj. 



B. Models with A x > g 



We now investigate what happens when A x 3> g. As discussed in the introduction, this occurs naturally in many 
models. For the run to be presented we used m 2 — 10~ 12 $q, g = 10~ 8 , A^ = 0, A x = 0.01, and Ax = 5n/2. The 
evolution in this case is very different from the A x = case (case 1(a)) considered above. As can be seen from figure 
6, the resonance starts growing just as in case 1(a). However, when the occupation numbers reach about 10 2 /A X , 
scatterings become fast and populate the infrared states. This gives rise to the backreaction Sm 2 , = 3A x (x 2 ), which 
shifts the resonances toward the infrared. As the first resonance sweeps through the infrared, it further populates 
the states. Once the resonance reaches k ~ 0, it is no longer very effective so that scatterings remove particles from 
the infrared faster than they are supplied by the resonance. This leads to depletion of the infrared states, which 
decreases the backreaction a little, increasing the resonance's effectiveness. This feedback mechanism leads to a state 
with slowly varying occupation numbers, just as in case 1(a). As a consequence, energy is removed from the zero 
mode at a slowly varying rate, starting at about t = 2.5 x 10 8 <I>q . In contrast to case 1(a) the <j> field plays no 
part in keeping the occupation numbers constant: it is the much faster x~ X scatterings which remove particles 
from the resonant momenta as they are created |0|. Consequently the EV does not get depleted by scatterings, and 



the rate of energy density loss of the </>-field is simply 



-p^ cs 2/ifeW0, which leads to an estimate of the decay 



time Tdecay 



lOw0<&5/(2/ifc/cjr es w rcs n£ ). The parameters can be deduced from our simulation. The 



relevant resonance in this case is the second one (recall that the first one has been shifted to the infrared) which 
can be seen at k « 0.07^/g^o at late times in figure 6. We estimate k 2 os ~ 0.5y/guj,p^o, /Zfc rcB ~ 0.03, the resonance 
width w Tes — 0.05fc rcs , and n£ ~ 15/A X , so that Tdecay ~ 10 4 /(ff n fc w<j>) ~ 10 15 $^ . Note that this time scale is 
parametrically larger than the resonance decay time by X x /glu.(l/g), which is about 10 5 in our case. We point out 
that (x 2 ) reaches a maximum value of about 6 x 10" 9 $g at t ~ 2 X lO 8 *^ 1 , and by t = 8 x lO 8 ®' 1 less than 0.0005% 
of the energy is in the x field fl9| |. The rate of energy density increase is about 5p x /St ~ 4 x 10~ 28 <&q. A linear 
extrapolation gives for the decay time Tdecay ~ P<p/(^Px/^) ~ lO 16 ^^" 1 , in rough agreement with our simple model. 
We emphasize that we do not claim that at very late times the simple feedback mechanism gives a correct picture. 
The main point of this simulation was to establish that a moderate self-coupling of the x field slows down the EV 
decay by many orders of magnitude. This conclusion remains valid when x couples moderately to other fields. The 
only difference is that x now decays into these fields, which in turn induces backreaction on x- While the details are 
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\<f> = 0, g = 10~ 8 and t* = y/g&ot). 



model dependent, the basic mechanism remains valid, and parametric resonance is rendered ineffective. 



IV. MASSLESS INFLATON (TYPE II MODELS) 

Let us now turn to type II models. These were investigated using the same techniques as above, so rather than 
present the results in detail we will simply state our conclusions. Before doing so it is worth pointing out that for 
these models the expanding universe equations of motion are conformally equivalent to those in Minkowski space. 
By this we mean that in terms of the variables r = J dt/a(t), <j> = (f>a(r) / a(0) , and \ — ^a(r)/a(0), where a(r) is 
the scale factor, the Friedmann-Robertson- Walker (FRW) equations of motion are of exactly the same form as those 
in ordinary static spacetime [po[ . Hence, with the simple replacements above, our numerical calculations for type II 
models apply directly to an expanding universe. 

The main results of our study are as follows. 



A. Models with A x = 

If X x — 0, the inflaton decays via parametric resonance in much the same way as it did in case 1(a). That is, after 
a brief period of exponential decay a slowly varying state sets in during which the decay is dominated by scatterings. 
For example, figure 7 shows the variances of the fields for m? — 0, A^ = 10~ 12 , and g — 1CP 8 . Since in this case 
ujcj, — Q.'&5^J~^,<&Q^ q w 3500. The figure indicates that, as in case 1(a), the exponential regime ends at t ~ 0.8 x 10 8 $q 1 . 
The maximum values of the variances are again to a good approximation given by g(x 2 ) ~ ^\ and g(4> 2 ) ~ u^^/g^o. 

The decay time can still be estimated from (Q), with the replacement — > ~ O.SSVA^o- The slowly varying 
value of gn£ is again about hto^j ^/g^o, so type I and II models behave very much alike for equal values of the 
parameters. One difference is that for type II models the inflaton energy decays somewhat faster, since ^ is larger 
and in addition there are scatterings via the X^cj) 4 term. For the A x = runs described in this work, the energy decays 
about 50% faster for the quartic potential. The value of Hk is about 25% larger for that case. 

As discussed at the beginning of this section, the lattice results for type II models can be mapped onto the 
expanding universe by the appropriate rescaling. To see how this works in practice consider the problem of obtaining 
the maximum value of (% 2 ) reached during preheating in the expanding universe. This value is interesting because 
it quantifies the effective temperature of non-thermal phase transitions Km . Since for type II models the universe is 
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FIG. 7. The variances of the fields as a function of time for a quartic inflaton potential (m = 0, A^, = 10 -12 , g = 10" 



radiation dominated we can write a(r)/a(0) = 1 + £/oo(0)t, where r is the conformal time and Hq = \f\ c j ) <& 2 / ^/Y2Mp 



is the Hubble constant at the end of inflation. As discussed at the beginning of section IV , in a Friedmann-Robertson- 
Walker universe our lattice fields are rescaled by a factor a(r)/a(0) and our lattice time is r. Hence (x 2 frw) = 
(a(0)/a(r)) 2 (x 2 ). Consider, for example, the variances shown in figure 7. Since (x 2 ) rises rapidly and then becomes 
slowly varying, (x 2 FR\v)max occurs at the moment when the fast growth terminates. This occurs at r w SOOO/y^'&o = 
8x lO 7 ^ 1 when ( X 2 ) ^6x10^. Hence we obtain (x 2 FRW ) max ~ 6 x 10~ 5 *§/(l + 2 x 1OVA $ O /Mp) 2 ~ 10" 7 M 2 , 
where we have used a(0) = 1 in our units. We can obtain general formulae for the maximum variances of the fields 
when q 3> 1 by recalling that the slowly varying scattering regime sets in when g((S4>) 2 ) ~ k 2 cs ~ y/gu$<5>o/2 and 
#(x 2 ) ~ where ~ O.85"\/A0<I>o. This occurs when n£ ~ {g^/q)~ x , i.e. at the time r ~ \n(l/g^/q)/2fik^rt>- 

Combining this with our equation for a(r) we obtain 

( W frw) 2 )™x ~ 0-2 1/21 2 * (Ha) 
q v ' A m (A^ q *t*) 

i 

—) 2 -3-^ Ml (lib) 

\gj ln 2 (A / 5 3 ) 

and 



(XfbwU ~ 0-06 : 2 , x -i 3/2 , m p (12a) 
gin (A g -V 2 ) 

— 1 -5-^ Af 2 , (12b) 

9 J ln 2 (A / ff 3) p 

where Mp as 2.4 x 10 18 GeV. These results can be considered upper limits for the variances since they were obtained 
for a massless \ field and a mass term in the equations of motion only suppresses the resonance Q. Note that 
(x 2 FRw)™aa; ^ 9 1 ( U P to logarithms), while Hartree type approximations which neglect scattering predict a q^ 1 / 2 
dependence 0. That the Hartree approximation may not give the correct q dependence of (% 2 ) was previously 
mentioned in 
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B. Models with A x > g 



If A x 3> g, the situation is drastically changed. Just as for type I models the decay into the x~ ne ld is extremely 
slow. But for type II models the inflaton can decay into its own fluctuations via parametric resonance, as discussed 
below Eq. (|J). This process, which is much slower than decay into x fluctuations for A x — 0, becomes dominant 
for A x ^> g. In fact the inflaton decays into its own fluctuations as if it were not coupled to the x _ne ld at all. We 
have verified this by running the one field case and comparing the results to the two field simulation with the same 
parameters in the <j> sector. For example, with m = 0, A^ = 10~ 12 , A x = 10~ 2 and g — 10~ 8 we find that 73% of the 
energy has decayed by t = 8 x 10 8 <f> " 1 and ((5</>) 2 ) max ~ 8 x 10 — 2< E>g, which agrees very well with the corresponding 
one field run. We also note that in the two field run (x 2 ) grows exponentially at first but then reaches a maximum 
value of about 3 x 10~ 9 <f )2 ) at t ~ 0.4 x lO 8 ^ 1 . The details of the one field case are discussed in ||. One may 
wonder why the perturbative energy flow via scatterings from <j) to x fluctuations is ineffective even after the infrared 
4> occupation numbers reach a huge value of order A^ 1 3> <7 _1 . The reason is that due to the backreaction g(S(f) 2 ) , x 
fluctuations become massive and <f>4> — > XX scatterings are kinamatically forbidden. (Indeed, taking the above value 
for ((£</>) 2 )max, one finds m x /uj <j , — 30.) 

We point out that it is very difficult to simulate type II models with large q on the lattice. The reason is that the 
typical resonant momenta for the two fields differ by a factor kf es /k* cs ~ 9 -1 ^ 4 , making it difficult to run simulations 
that capture both resonances with sufficient infrared and ultraviolet resolution. Our numerical results on 128 3 lattices 
were hence obtained in a two step process: we first did preliminary runs, choosing Air optimally for each resonance in 
turn. This allowed us to determine which one was dominant, namely the x resonance for A x = and the <fi resonance 
for X x g. We then did extended runs choosing Ax just large enough to capture the physics of the dominant 
resonance, giving us the best possible ultraviolet range for the given situation. 



V. CONCLUSION 



We have numerically investigated the decay of the inflaton coupled to a massless scalar field. Our main results for 
both quadratic (type I) and quartic (type II) inflaton potentials are as follows: 

1. If the self-coupling of the decay product is small (A x <C g) we find two distinct stages of inflaton decay for q > 25: 
an exponential regime in which the field decays via parametric resonance, followed by a scattering dominated regime. 
The fraction of energy which decays by the time the scattering regime begins is of order q~ x / 2 . This means that for 
5> 1 scatterings are responsible for most of the inflaton decay, and we find that the decay time scale in the scattering 
regime is significantly longer than the resonance decay time scale. 

2. If X x 3> g (which is natural for many models) we find that parametric resonance cannot transfer energy to x- The 
reason is that scatterings due to the self-interactions limit the maximum occupation numbers of x, allowing only a tiny 
fraction of the inflaton energy to be transferred during the resonant stage. In the subsequent scattering dominated 
regime the transfer of energy is extremely slow, leading to huge decay times for type I models. For type II models 
the inflaton decays into its own fluctuations, essentially as if it were not coupled to the second field at all. Although 
these results were obtained in a simple model with two scalar fields, the basic mechanism, and hence our conclusion, 
remains valid for a realistic model with many mutually interacting fields. 

In this work we have studied inflaton decay by integrating the Minkowski space equations of motion. It is an 
interesting question how our conclusions might change in an expanding universe. For type II models the answer is 
simple: as explained in secti on frv| , our analysis is directly applicable after the appropriate rescaling. The mapping 



was made explicit in section IV A, where we obtained upper limits for the variances of the fields in the expanding 
universe. 

For type I models the situation is more complex. The time scale H^ 1 defined by the the expansion rate H 



VAa/3/Mp ^5 ^</./2 is much shorter than the scattering decay time scale (see Eq. (10)). Consequently we cannot 
expect our work for type I models to be a good approximation to the expanding universe case. However, we do expect 
our conclusion that the inflaton cannot decay into a field with moderate self-interactions to remain valid. First of 
all, parametric resonance in an expanding universe is suppressed. Second, in the scattering regime the maximum 
occupation number will remain of order A" 1 , so one expects that the total energy deposited in the x field fluctuations 
remains a tiny fraction of the inflaton energy even in an expanding universe. 
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